package fem2d;

import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IVector;

import java.util.ArrayList;
import java.util.Collections;

import math2.ConstantFunctionRnToR;
import math2.FunctionRnToR;
import fem2.AbstractStructuralStaticDemo;
import fem2.Constraint;
import fem2.DistributedLoad;
import fem2.Edge;
import fem2.Element;
import fem2.Load;
import fem2.MaterialModel;
import fem2.Mesh;
import fem2.MeshGenerator;
import fem2.MeshPart;
import fem2.MeshUtilities;
import fem2.Model;
import fem2.Node;
import fem2.Observer;
import fem2.StringUtilities;
import fem2.element.StructuralElement;
import fem2.enu.SortingCriteria;
import fem2.enu.State;
import fem2.enu.StressType;
import fem2.material.StVenantKirchhoffMM;
import fem2.material.SPRrecoveryMM;
import fem2.observer.GidStressObserver;
import fem2.observer.TimeStepObserver;
import fem2.pre_and_post.GidMeshGenerator;
import fem2.strategies.Strategy;

public class PerforatedPlate extends AbstractStructuralStaticDemo {

	public PerforatedPlate(String meshFileName) {
		projectDir = "/home/hbui/kratos_janosch";
		projectName = "perforated_plate.gid";
		this.meshFileName = meshFileName;
	}

	@Override
	public Mesh createMesh() {
		String fn = projectDir + '/' + projectName + '/' + meshFileName;
		MeshGenerator mg = new GidMeshGenerator(fn);
		return mg.getMesh(2);
	}

	@Override
	public Model createConditions(Model m) {
		Mesh mesh = m.getMesh();

		ArrayList<Node> bottom = MeshUtilities.seekNodesOnSurface(mesh, 0, 1, 0);
		ArrayList<Node> left = MeshUtilities.seekNodesOnSurface(mesh, 1, 0, 0);
		ArrayList<Node> right = MeshUtilities.seekNodesOnSurface(mesh, 1, 0, -4);

		Constraint cx = new Constraint(false, true);
		Constraint cy = new Constraint(true, false);

		mesh.addConstraint(cy, bottom);
		mesh.addConstraint(cx, left);

		mesh.setSortingCriteria(SortingCriteria.SORT_ON_Y_COORDINATE, right);
		Collections.sort(right);
		mesh.setSortingCriteria(SortingCriteria.SORT_ON_ID, right);

		ArrayList<Edge> right_edges = MeshUtilities.createEdges(right);
		FunctionRnToR F = new ConstantFunctionRnToR(1.0);
		for (Edge e : right_edges) {
			mesh.addEdge(e);
			Load l = new DistributedLoad(e, F, null);
			m.addLoad(l);
		}

		return m;
	}

	private IVector computeExactStress(double x, double y) {
		double r = Math.sqrt(Math.pow(x, 2) + Math.pow(y, 2));
		double theta = Math.atan2(y, x);
		double a = 1.0;
		double a2 = Math.pow(a, 2);
		double a4 = Math.pow(a, 4);
		double r2 = Math.pow(r, 2);
		double r4 = Math.pow(r, 4);
		double cos2 = Math.cos(2 * theta);
		double cos4 = Math.cos(4 * theta);
		double sin2 = Math.sin(2 * theta);
		double sin4 = Math.sin(4 * theta);

		double oxx = 1 - (a2 / r2) * (1.5 * cos2 + cos4) + 1.5 * (a4 / r4) * cos4;
		double oyy = -(a2 / r2) * (0.5 * cos2 - cos4) - 1.5 * (a4 / r4) * cos4;
		double txy = -(a2 / r2) * (0.5 * sin2 + sin4) + 1.5 * (a4 / r4) * sin4;

		return new ArrayVector(oxx, oyy, txy);
	}

	@Override
	public MaterialModel createMaterial(Model m) {
		double E = 1e3;
		double nu = 0.3;
		double t = 1.0;

		MaterialModel mm = new StVenantKirchhoffMM(E, nu, t, 0, State.PLANE_STRAIN) {

			@Override
			public IVector getStress(Element e, int k, StressType type) {
				switch (type) {
				case EXACT_STRESS:
					double[] xi = e.getIntegrationPoint(k);
					double[] X = e.getMeshPart().XAt(xi);
					return computeExactStress(X[0], X[1]);
				case STRESS_ERROR:
					xi = e.getIntegrationPoint(k);
					X = e.getMeshPart().XAt(xi);
					IVector error = computeExactStress(X[0], X[1]);
					IVector s = ((StructuralElement) e).getStress(k);
					BLAM.add(-1.0, s, error);
					return error;
				default:
					return super.getStress(e, k, type);
				}
			}

		};

		return new SPRrecoveryMM(m, mm);
	}

	@Override
	public Element createElement(MeshPart mp, MaterialModel mm) {
		return new StructuralElement(mp, mm);
	}

	@Override
	public void addObservers(Model m, Strategy s, TimeStepObserver o) {
		super.addObservers(m, s, o);

		String meshName = StringUtilities.getFileName(meshFileName);

		Observer so1 = new GidStressObserver(m, o, projectDir, projectName, meshName,
				StressType.EXACT_STRESS);
		Observer so2 = new GidStressObserver(m, o, projectDir, projectName, meshName,
				StressType.STRESS_ERROR);
		Observer so3 = new GidStressObserver(m, o, projectDir, projectName, meshName,
				StressType.RECOVER_STRESS_ZZ);

		s.addObserver(so1);
		s.addObserver(so2);
		s.addObserver(so3);
	}

	public static void main(String[] args) {
		// AbstractStructuralStaticDemo demo1 = new
		// PerforatedPlate("T3_0.5.msh");
		// demo1.run();
		// AbstractStructuralStaticDemo demo2 = new
		// PerforatedPlate("T3_0.25.msh");
		// demo2.run();
		// AbstractStructuralStaticDemo demo3 = new
		// PerforatedPlate("T3_0.1.msh");
		// demo3.run();
		// AbstractStructuralStaticDemo demo4 = new
		// PerforatedPlate("T3_0.05.msh");
		// demo4.run();

		AbstractStructuralStaticDemo demo1 = new PerforatedPlate("Q4_0.5.msh");
		demo1.run();
		AbstractStructuralStaticDemo demo2 = new PerforatedPlate("Q4_0.25.msh");
		demo2.run();
		AbstractStructuralStaticDemo demo3 = new PerforatedPlate("Q4_0.1.msh");
		demo3.run();
		AbstractStructuralStaticDemo demo4 = new PerforatedPlate("Q4_0.05.msh");
		demo4.run();
	}
}

/**
 * results:
 * 
 * T3_0.5: <br>
 * Show Min Max 'Sxx-STRESS_ERROR_TENSOR': Min = -0.83826, Max = 0.22607<br>
 * Show Min Max 'Syy-STRESS_ERROR_TENSOR': Min = -0.33456, Max = 0.43576<br>
 * Show Min Max 'Sxy-STRESS_ERROR_TENSOR': Min = -0.093075, Max = 0.11428<br>
 * 
 * T3_0.25<br>
 * Show Min Max 'Sxx-STRESS_ERROR_TENSOR': Min = -0.86171, Max = 0.27786<br>
 * Show Min Max 'Syy-STRESS_ERROR_TENSOR': Min = -0.40888, Max = 0.58531<br>
 * Show Min Max 'Sxy-STRESS_ERROR_TENSOR': Min = -0.10155, Max = 0.13945<br>
 * 
 * T3_0.1<br>
 * Show Min Max 'Sxx-STRESS_ERROR_TENSOR': Min = -0.79954, Max = 0.29375<br>
 * Show Min Max 'Syy-STRESS_ERROR_TENSOR': Min = -0.43375, Max = 0.56946<br>
 * Show Min Max 'Sxy-STRESS_ERROR_TENSOR': Min = -0.097947, Max = 0.16152<br>
 * 
 */
